  library(mlogit)
  library(plyr)
  library(reshape2)
  #library(lme4)
  library(rbenchmark)
  library(Matrix)
  library(MatrixModels)
  
  setwd("~/Research/2014/rmm")
  source("R/classes_methods.R")
  source("R/utilities.R")
  source("R/getFrame.R")
  source("R/getData2.R")
  dyn.load("src/rmm.so")
  .Call("init")


#'
#' The high dimensional logit model using sparse matrix
#'
#' @param formula the input model formula
#' @param data input data of class mlogit 
#' @param nrho number of rhos to specify the grid  
#' @param rho.range the min and max of rhos that are used to specify the grid
#' @param rho a numeric vector. If not null, this defines the values of rhos 
#' to be used in the estimation
#' @param max.nbr the maximum number of neighbors to search for
#' @param sigma the standard deviation in the Gaussian kernal 
#' @param print whether iteration info should be printed out
#'
#' @return  an object of class hlogit 
#'

rmm <- function(formula, data, model = c("nrmm", "elmm", "glmm"), 
                ncluster = 1, subset, weights = NULL, 
                nlambda = 20, lambda = NULL, 
                penalty.factor = NULL, 
                doFit = TRUE, inits = NULL, control = rmm.control()){  
  
  mc <- match.call()
  model <- match.arg(model)  
  if (ncluster > 1 && !is.null(weights))
    stop("'weights' are not supported for mixture models.")
  #check.data(data)
  # get design matrix
  fr <- getFrame(mc, formula, control["use.contrast"])
  X <- fr$X  
  if (model == "glmm")
    control["standardize"] <- TRUE
  if (control["standardize"]){
    sx <- scaleX(fr$X, fr$gvar)
    X <- sx$X
    attr(X, "U") <- sx$XU
  }    
  
  # get lambda
  intercept <- attr(terms(fr$mf), "intercept")
  if (model == "nrmm"){
    nlambda <- 1
    lambda <- 0
  } else {
    if (is.null(lambda)){    
      maxLam <- mnl.get.maxLam(X, fr$y, fr$gvar, intercept)    
      maxLam <- maxLam * control["lambda.max.scale"]
      lambda <- getLambda(nlambda, maxLam, control["lambda.min.ratio"])      
    } 
    nlambda <- length(lambda)
  }
  # the dims slot
  dims <- as.integer(c(nrow(X), ncol(X), nlevels(fr$chid), nlevels(fr$alt),  
                       length(fr$gvar) - 1, nlambda, max(fr$id) + 1, ncluster))
  names(dims) <- c("nO", "nB", "nP", "nC", "nG", "nL", "nI", "nS")  
  
  # get the penalty factor
  pfct <- penalty.factor
  if (is.null(pfct)){    
    if (intercept){
      pfct <- c(0, rep(1, dims["nG"] - 1))
    } else {
      pfct <- rep(1, dims["nG"])
    }
  }
  
  L <- new("dCHMsimpl")      
  if (model == "nrmm" && !control["newton"]){
    xtx <- crossprod(X)
    L <- Cholesky(xtx, super = FALSE, LDL = FALSE)
  }
  struct <- new("mm", call = mc, X = X, y = fr$y,
                dims = dims, wts = fr$wts,
                control = control, frame= fr$mf,
                subject = fr$id, gid = fr$gid,                
                gvar = fr$gvar, formula = formula)
  ans <- switch(model, 
                nrmm = new("nrmm", struct, L = L),
                glmm = new("glmm", struct, 
                           lambda = as.numeric(lambda),
                           pfct = as.numeric(pfct)))
  if (!is.null(inits)){
    if (length(inits) != dims["nB"] * dims["nS"])
      stop("wrong length of inits")
    ans@beta[, 1] <- inits
  }
  # return the object without fitting 
  if (!doFit) return(ans)        
  
  # optimize 
  switch(model, 
         nrmm = .Call("R_do_nrmm", ans),
         glmm = .Call("R_do_glmm", ans),
         elmm = .Call("R_do_elmm", ans))
  
  if (ans@control["cvg"] != 0)
    warning("Algorithm does not converge!")
  
  return(ans)  
}

